program flu2d 
  use technical
  use mod_2dflu
  use mod_part_interp
  use mod_read_write
  use immersed_boundary
  use force
  use fft
  use omplib
  implicit none
!! other defination of variables 
  integer :: i1,i2,i3,ispectra,iouter,isnap,nouter,icall,st
  double precision :: En,Ent,Entdiss,epsilon,sumEkbyki,compare,dummy,tstart,tend
  integer :: k1,k2,ksqr,ireal,iimag,mshl,iloop
  logical :: lfft_in_place=.true.
  integer :: Np=0,nstokes=0,nthreads
  integer :: ierr
  double precision :: time
  double precision :: allocated_memory=0.
  character (len=labellen) :: fname1,fname2
!  character*500 :: cnpar,formp
!! --------------------------------------------------------------
  write(*,*) '*****************************************************************'
  write(*,*) '     SPECTRAL SOAP   '
  write(*,*)'2-d Spectral Code for Hydrodynamics (Incompressible Navier-Stokes)'
  write(*,*) 'Reading input parameters ...'
!! % first read the basic variables
  if(file_exists('para_2d.in')) then 
    open(unit=1,file='para_2d.in',status='old')
    call read_basic(1,ierr)
    if (ierr.ne.0) call fatal_error('flu2d','cannot read namelist basic_param')
    if (lforcing)  then 
      call read_force_pars(1,ierr)
      if (ierr.ne.0) call fatal_error('flu2d','cannot read namelist force_param')
    else
      lfourier_forcing=.false.;lrspace_forcing=.false.
    endif
    if (lparticle) call read_particles_param(1,ierr)
    if (ierr.ne.0) call fatal_error('flu2d','cannot read namelist particle_param')
    if (limmersed_boundary) call read_immersed_boundary(1,ierr)
    if (ierr.ne.0) call fatal_error('flu2d','cannot read namelist immersed_bound_param')
  else
    call fatal_error('input file para_2d.in','not found')
  endif
 write(*,*) '..... done'
 write(*,*) 'Resolution:',nn
 write(*,*) 'viscosity:',vis
 write(*,*) 'Fixed time step=',delta
 write(*,*) '-------------------------------------------------'
!! ---------------------------------------------!!
! check if data directory exist. If not, then stop. 
 if( .not.file_exists('data')) call fatal_error('main:','datadir (typically data) directory must exist')
!$OMP PARALLEL 
!$OMP MASTER
 call get_total_threads(nthreads) 
 if (nthreads .gt. 1) then
    write(*,*) 'Using threads with number of threads = ',nthreads
 else
    write(*,*) 'Not using any threads, nthreads = ',nthreads
 endif
!$OMP END MASTER
!$OMP END PARALLEL
 write(*,*) '-------------------------------------------------'
 if (lcalc_Gradu) write(*,*) 'Gradu Matrix will be calculated '
 write(*,*) 'Initializing fluid .......'
 !! %---------------------
 !! \section{Initialisation}  
 call intialize_grid()
 write(*,*) '............. done'
 write(*,*) 'Now allocating fluid arrays......... '
!! allocate and calculate the global arrays --------
  allocate(den_state(0:nshell))
  allocated_memory=allocated_memory+4.*dfloat(nshell+1)
  allocate(time_increment(0:ksqr_max),alpha(0:ksqr_max))
  allocated_memory=allocated_memory+2.*8.*dfloat(ksqr_max+1)
  den_state = 0
  time_increment = 0.0d0
  call global2d
!! %---and also energy and dissipation arrays ------
  allocate(E_Omega(0:nshell,3))
  E_Omega = 0.0d0
  allocated_memory=allocated_memory+3.*8.*dfloat(nshell+1)
!! %-stream-function and vorticity array are allocated  ----
  allocate(psi(n1+2,n2),omega(n1+2,n2))
  allocated_memory=allocated_memory+2.*8.*dfloat((n1+2)*n2)
  if(lcalc_Gradu) then
    allocate(dx_ux(n1+2,n2),dy_ux(n1+2,n2),dx_uy(n1+2,n2))
    allocated_memory=allocated_memory+3.*8.*dfloat((n1+2)*n2)
    dx_ux=0.;dy_ux=0.;dx_uy=0;
  endif
!! %-local arrays to do transformations  ------------------
  psi = 0.0d0
  omega = 0.0d0
!! %-allocate Jacobean etc 
  allocate(jac_old(n1+2,n2))
  allocated_memory=allocated_memory+8.*dfloat((n1+2)*n2)
  jac_old = 0.0d0
!! -- allocate velocity in fourier space 
  allocate(ukx(n1+2,n2),uky(n1+2,n2)) 
  allocated_memory=allocated_memory+2.*8.*dfloat((n1+2)*n2)
  ukx=0.;uky=0.
  write(*,*)'..........done'
  write(*,*)'Allocated memory (in GB) so far is:',allocated_memory/(10.**9)
 !! --- Allocate arrays for particle tracking
 if (Np.eq.0) then
   lparticle=.false.
   write(*,*) 'ZERO particles: No particles allocated'
 else
   lparticle=.true. 
   write(*,*)'Now allocating memory for particles ....'
   call allocate_particles(Np,nstokes,lcalc_Gradu,allocated_memory)
 endif
 write(*,*)'..........done'
 if (limmersed_boundary) then
   call allocate_immersed_boundary()
 else
   write(*,*) 'Not using Immersed Boundaries'
 endif
 write(*,*)'Allocated memory (in GB) so far is:',allocated_memory/(10.**9)
  !
  !! -------------------------------------------------
  !! Initialization for Fourier Transform
  !! -------------------------------------------------
  write(*,*) 'Initializing Fourier Transform ...'
  dim(1) = n1 
  dim(2) = n2  
  scale = 1.0d0/dfloat(n1*n2)
  lfft_in_place=.true.
  call initialize_2dfft(dim,lfft_in_place)
  write(*,*) ' ........done'
!!
!! %-------------------******************------------------------------
!! we set an initial configuration.  
  if(nrun.eq.1) then
    write(*,*) 'Setting up Initial condition for velocity....'
    call init_uu(oamp,kini)
! ---------------
  else
!! Otherwise (i.e. if this is not the first run ) we read the
!! configuration of vorticity
    write(*,*)'Reading real space vorticity data from previous run...'
    if(file_exists('vortex.in')) then 
      open(unit=11,file='vortex.in',status='old')
      do i2=1,n2
        read(11,*) (omega(i1,i2),i1=1,n1)
      enddo
      close(11)
    else
      call fatal_error('nrun .ne. 1','but input file vortex.in not found')
    endif
!    open(unit=11,file='input/test/vor_real2163.dat',form='unformatted',status='old')
!     open(unit=11,file='input/vortex.out',form='unformatted',status='old')
!    read(11) ((omega(i1,i2),i1=1,n1),i2=1,n2)
!    close(11)
!     open(unit=11,file='input/vor.2163',form='unformatted',status='old')
!     read(11) ((omega(i1,i2),i1=1,n1+2),i2=1,n2)
!     close(11)
!     do i1=n1+1,n1+2
!        omega(i1,:) = 0.0d0
!     enddo
  endif
  write(*,*) '.....done'
! Now write the initial condition to a file
  call write_2dScalarField_ascii(omega(1:n1,1:n2),'oini',n1,n2)
! 
!After either reading in vorticity or generating it from intial configuration
!Calculate velocity, spectra, etc
!But, first vorticity to Fourier space
!
  write(*,*)'FFT vorticity to Fourier space...'
  call fft_2d_real2complex(omega)
  write(*,*) '...........done'
!! obtain stream function and velocity from vorticity
  write(*,*)'Calculating velocity, psi and energy from vorticity ...' 
  call getpsivel
  call energy2d
  write(*,*)'........done, now writing initial energy...'
!
! write them out.
!
    open(unit=110,file='initial_spectrum.out',status='unknown')
    do ispectra = 1,nshell
      write(110,*) dfloat(ispectra),E_Omega(ispectra,1)
    enddo
    close(110)
  write(*,*)'............done'
!  ------------- Initialize the particle positions
! But initializing particles might need ux and uy in real space
!
    write(*,*)'Calculating velocity in real space (two FFTs) ...'  
    call fft_2d_complex2real(ukx)
    call fft_2d_complex2real(uky)
    do i1=n1+1,n1+2
      ukx(i1,:) = 0.0d0
      uky(i1,:) = 0.0d0
    enddo
    ukx = ukx*scale
    uky = uky*scale
    ux=>ukx(1:n1,1:n2)
    uy=>uky(1:n1,1:n2)
    write(*,*)'..........done'
!-------------------
    if (nprun.eq.1) then
      write(*,*)'Initializing particles ........'
      call init_particle(n1,n2,dx,dy,ux,uy)
    else     
      write(*,*)'Reading particle data from previous run ....' 
      call read_input_particle()
    endif
    write(*,*) '.........done'
!-------------------
    if (limmersed_boundary) then
      write(*,*) 'Initializing immersed boundary ...'
      call initialize_immersed_boundary()
      write(*,*) '... Done'
    endif
!-------------------
    if (lforcing) then
      write(*,*) 'Initializing external force ...'
      call initialize_force()
      write(*,*) '... Done'
    endif
!-------------------
    write(*,*) '****************************************************'
    write(*,*) 'Now starting time marching:     '
!! ----------------------------
  !! \section{Time Marching}
  open(unit=15,file='energies.out',status='unknown')
!  open(unit=17,file='balance.out',status='unknown')
  call open_output_particle()
!!
  call wall_clock_time(tstart)
  time=0.
  icall=0
!  do icall = 1,maxiter/navg
   do while ( (icall .lt. maxiter/navg) .and. (time.lt.tmax) )
    do iinner = 1,navg
      do iloop = 1,nloop
         call rnkt(mod(iloop,nloop))
!         write(*,*)'inner loop,middle loop, outer loop =',iloop,iinner,icall
         time=time+delta
         icall=icall+1
       enddo
       if( file_exists('STOP',.true.)) call stop_it('found STOP file')
       call getpsivel
       E_Omega = 0.0d0
       call energy2d
       write(15,*) 0.5d0*sum(E_Omega(:,1)),0.5*sum(E_Omega(:,2))
       write(*,*) 'time,Energy,Enstrophy=',time,0.5d0*sum(E_Omega(:,1)),0.5*sum(E_Omega(:,2))
!      write(17,*) fw,fu,sum(E_Omega(:,3))
!     write(fname2,'("data/lambda.",I4.4)')(iinner+(icall-1)*navg)
!     open(unit=333,file=TRIM(fname2),form='unformatted',status='unknown')
!     write(333)((lambda(i1,i2),i1=1,n1,4),i2=1,n2,4)
!     close(333)
!     write(fname2,'("data/psi.",I4.4)')(iinner+(icall-1)*navg)
!     open(unit=333,file=TRIM(fname2),form='unformatted',status='unknown')
!     write(333)((psi(i1,i2),i1=1,n1,4),i2=1,n2,4)
!     close(333)
	!
       if(write_fullparticle)then
         call write_output_particle
       else
         call write_lambda_particle
       endif
     enddo
     write(fname2,'("spectras/spectra.",I4.4)')(icall)
     open(unit=110,file=TRIM(fname2),status='unknown')
     do ispectra = 1,nshell
       write(110,*) dfloat(ispectra),E_Omega(ispectra,1),E_Omega(ispectra,2)
     enddo
     close(110)
!! % ----Write the streamfunction and the vortictiy
!     write(fname2,'("data_psi/psi.",I4.4)')(icall)
!     open(unit=333,file=TRIM(fname2),form='unformatted',status='unknown')
!     write(333)((psi(i1,i2),i1=1,n1+2),i2=1,n2)
!     close(333)
     call write_trajectory_particle
     !!
   enddo
   call wall_clock_time(tend)
   write(*,*) 'The code has ran for total',maxiter*nloop ,'timesteps'
   write(*,*) 'In time (in sec)', (tend-tstart)
   write(*,*) 'with number of threads =',nthreads
   write(*,*) 'Time taken(per grid point, per time step) = ', (tend-tstart)/(dfloat(n1*n2)*(maxiter*nloop))
   close(15)
   ! close(17)
   call close_output_particle
!!End of time marching
!! ------------------- write output ------------------------- !!
   if (lwrite_velocity_output) call getpsivel
! this gives velocity in fourier space, if we choose to write it in
! real space we need to fourier transform:
   if (lwrite_vel_rspace) then
     call fft_2d_complex2real(ukx)
     call fft_2d_complex2real(uky)
     ukx = ukx*scale
     uky = uky*scale
     do i1=n1+1,n1+2
       ukx(i1,:) = 0.0d0
       uky(i1,:) = 0.0d0
     enddo
     if (lascii_output) then
       call write_2dScalarField_ascii(ukx(1:n1,1:n2),'ux',n1,n2)
       call write_2dScalarField_ascii(uky(1:n1,1:n2),'uy',n1,n2)
     else
       call write_2dScalarField_binary(ukx(1:n1,1:n2),'ux',n1,n2)
       call write_2dScalarField_binary(uky(1:n1,1:n2),'uy',n1,n2)
     endif
   else
     if (lascii_output) then
       call write_2dScalarField_ascii(ukx,'ukx',n1+2,n2)
       call write_2dScalarField_ascii(uky,'uky',n1+2,n2)
     else
       call write_2dScalarField_binary(ukx,'ukx',n1+2,n2)
       call write_2dScalarField_binary(uky,'uky',n1+2,n2)
     endif
   endif
!! Save vorticity file
   call fft_2d_complex2real(omega)
   omega = omega*scale
   do i1=n1+1,n1+2
     omega(i1,:) = 0.0d0
   enddo
   if (lascii_output) then
     call write_2dScalarField_ascii(omega(1:n1,1:n2),'omega',n1,n2)
   else
     call write_2dScalarField_binary(omega(1:n1,1:n2),'omega',n1,n2)
   endif
!! Also save the streamfunction if needed
   if (lwrite_streamfunction) then
     call fft_2d_complex2real(psi)
     psi = psi*scale
     do i1=n1+1,n1+2
       psi(i1,:) = 0.0d0
     enddo
     if (lascii_output) then
       call write_2dScalarField_ascii(psi(1:n1,1:n2),'stream',n1,n2)
     else
       call write_2dScalarField_binary(psi,'stream',n1,n2)
     endif
   endif
!!-----------------------------------
!! Finally we deallocate the arrays allocated,  
   deallocate(den_state)
   deallocate(time_increment)
   deallocate(psi,omega,ukx,uky)
   deallocate(E_Omega)
!! Destroy the plans created 
   call fft_2d_cleanup()
   write(*,*) 'End'
  !! %-------------*********************-------------------------------
  !! and end our program 			
 end program  flu2d
